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Abstract. We present the results of an exhaustive study of the ionized gas in NGC 588, a giant Hn region in 
the nearby spiral galaxy M33. This analysis uses a high number of diagnostics in the optical and infrared ranges. 
Four temperature diagnostics obtained with optical lines agree with a gas temperature of 11000K, while the 
[Om] A5007/A88^im ratio yields a much lower temperature of «8000K. This discrepancy suggests the presence of 
large temperature inhomogeneities in the nebula. We investigated the cause of this discrepancy by constructing 
photoionization models of increasing complexity. In particular, we used the constraints from the Ha and H/3 
surface brightness distributions and state-of-the-art models of the stellar ionizing spectrum. None of the successive 
attempts was able to reproduce the discrepancy between the temperature diagnostics, so the thermal balance of 
NGC 588 remains unexplained. We give an estimate of the effect of this failure on the O/H and Ne/O estimates 
and show that O/H is known to within ±0.2 dex. 
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1. Introduction 

Giant Hn regions (GHRs) are the most popular tracers 
of elemental abundances in distant galaxies, since they 
are easy to observe and, in principle, easy to analyze. 
Abundances can be obtained directly from the intensities 
of emission lines, without having to go through a detailed 
mode l analysis to extract the information fe.g.. lStasinskal 
However, for metal-rich Hn regions - say, with an 
oxygen abundance larger than half that of the Sun - the 
methods for abundance derivation are statistical and rely 
on previous calibrations based on photoionization mod- 
els. The large number of calibrations that have been pub- 
lished, leading to sign ificantly different abundances (e.g., 
iKennicutt et alJ |2003), show that the interpretation of 
emission lines in Hn regions may not be as easy as one 
would like. In the case of metallicities significantly below 
solar, abundance determination is considered more reli- 
able, since the optical emission lines provide a direct mea- 
sure of the gas temperature, which allows one to com- 
pute the elemental abundances directl y from the obser ved 
line intensities. However, as shown bv iPeimberti l|1967(l . if 
the temperature in Hn regions is not uniform but instead 
presents spatial fluctuations, the derived abundances will 
be biased. The reality of the presence of these temperature 



fluct uations and th eir cause are still a matter of debate 
(e.g.. lEsteba3l2002|) . 

Comprehensive analysis of nearby giant Hn regions are 
necessary to check our understanding of the thermal struc- 
ture of these objects and to validate empirical methods for 
abundance determinations. There have already been de- 
tailed photoionization stud i es of some giant Hii regions 
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most cases, the models were not able to reproduce the ob- 
served temperature indicators correctly. However, this was 
not the main goal of most of those studies, and the degree 
of sophistication of the models was perhaps insufficient for 
that purpose. 

In the present paper, we propose a comprehensive anal- 
ysis of NGC 588, a giant Hn region located on the out- 
skirts of the nearby spiral galaxy M33, with the aim of 
understanding its temperature structure. We gathered a 
large set of spectroscopic and imaging data in various 
wavelength ranges. From this set of data, we were able 
to give a full description of the ioni zing stellar popula - 
tion, based on a star- by-star analysis l|jamet et al.ll200^ . 
We now use the ionizing radiation field from this popula- 
tion together with information on the nebular morphology 
given by narrow band images to construct photoioniza- 
tion models. Constraints are provided by the strengths of 
optical and infrared emission lines. The availability of in- 
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Table 1. Journal of Observations. 



Aa (arcsec) 

Fig. 1. Ha image with the spectroscopic slits overplotted. 
North is up, and east is on the left. The image display 
scale is logarithmic. 



frared data is particularly important, since they enlarge 
the number of possible spectral diagnostics. 

The paper is organized as follows. In Sect. |21 we de- 
scribe the observational data and their processing. In Sect. 
|3]we present empirical diagnostics of the density tempera- 
ture, and chemical composition. In Sect. 0] we give details 
on the model-fitting strategy that we adopted. Starting 
from very simple models (Sect.EJ), we gradually increased 
the degree of sophistication in order to match the observa- 
tions as closely as possible (Sect. EH7J). We then examine 
the effects of energy sources other than the ionizing flux 
of the cluster (Sect. [HJ) . In Sect. OH we evaluate the impact 
of the unknowns of the thermal structure of the nebula 
on the determination of the O/H and Ne/O abundance 
ratios. Finally, we present our conclusions in Sect. 1101 

2. Observational data and reduction 

We acquired ground-based optical spectra of NGC 588 and 
retrieved a series of data available for this object in the 
archives of the Infrared Space Observatory (ISO) and of 
the Isaac Newton Group (ING). The main properties of 
the data with which we realized measurements are sum- 
marized in Table m Figure shows the nebula in Ha (see 
Sect. I2.2fl on which the spectroscopic slits were superim- 
posed. 

2.1. Optical CAHA spectra 

We obtained a set of 3 spectra on 27 August 1998 with 
the 3.5m telescope of the Centro Astronomico Hispano 
Aleman (CAHA, Calar Alto, Almeria, Spain). These ob- 
servations were acquired with the dual beam TWIN spec- 
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CAHA (PA=121°, width=1.2", 
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4 x 1800 
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2 x 1800 
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2 x 1800 


ISO/SWS (PA 


=-29°, data set 81601775, 1998-02-08) 


line 


range (/im) 


aperture (") 


[Siv] All/im 


10.45-10.55 


14x20 


[Nem] A16/mi 


15.46-15.64 


14x27 


[Sin] A19jum 


18.61-18.77 


14x27 


[Oiv] A26^m 


25.80-26.08 


14x27 


[Sin] A34/im 


33.17-33.75 


20x33 


[Nem] A36/im 
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Data set 


Spectrometer 


Exposure (s) 
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81601776 
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1402 


Imaging 
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Band 


Filter 


Exposure (s) 


H/3 


4861/54 


2 x 1800 


Ha 


6563/53 


1800 


Ha continuum 


6834/51 


1300 


[Om] A5007 


5007/51 


1800 



trograph, equipped with a beam splitter centered on 5500 
A which sends the signal through two arms ("red" and 
"blue"). The spectrum observed with each arm was col- 
lected with an SITe-CCD of 2000x800 15 /mi pixels. We 
took four exposures in the range 3595-5225 A (range B) 
with the blue arm, two in the range 5505-7695 A (range 
Rl), and two in the range 7605-9805 A ( range R2) with 
the red arm; each of these exposures was 1800 s long. The 
two gratings used for the observations, T7 (blue arm) and 
T4 (red arm), give dispersions of 0.81 and 1.08 A/pix, re- 
spectively. The signal was acquired through a 240 x 1.2" 
slit oriented at PA=121° with a resolution of ^2 A along 
the total spectral range, as measured from the FWHM of 
the lines of the sky background. 

2.1.1. Reduction 

We reduced the bidimensional spectra following the stan- 
dard procedure, mainly making use of the IRAF 1 package 
noao . twodspec . longslit. 

The flat-field correction was made with two light 
sources: an internal tungsten lamp for removal of the ef- 
fect of pixcl-to-pixel sensitivity variations and the twilight 
sky for correction for instrumental vignetting. Though the 
tungsten lamp spectra were heavily affected by fringing ef- 



1 IRAF is distributed by the National Optical Astronomy 
Observatory, which is operated by the Association of 
Universities for Research in Astronomy, Inc., under cooper- 
ative agreement with the National Science Foundation. 
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fects, the overal sensitivity flattening was satisfactory; we 
estimated the residuals to be ~1%. 

The wavelength calibration was performed from spec- 
tra of a helium-argon arc lamp, typicall givingy 20 useful 
lines per frame. Several of these spectra were obtained 
throughout the night in order to account for instrumental 
drifts. 

No data was available to calibrate the frames with re- 
spect to the positions along the slit, and the inspection 
of several spectra of stars showed ^2 pixel (« 1") dis- 
torsions around horizontal lines of pixels. Since we were 
mainly interested in integrated line fluxes, we ignored this 
issue, except for the hydrogen lines used for dereddening 
(see Sect. ETUI . 

We carefully combined the individual exposures in 
each range (B, Rl, R2). We first performed slight shifts 
of these frames along the slit axis, to match the profiles 
of the lines, since the instrument centering on the object 
was subject to small variations between the different ex- 
posures. The differences between the coregistered frames 
showed no detectable residuals other than cosmic rays, 
except for some of the brightest lines (Ha, H/3, and [Om] 
A5007) in a ~4" zone around the maximum nebular emis- 
sion. However, these discrepancies were small and resulted 
in flux errors less than ~3% in this zone. Since this error 
is small and concerns only a small zone of the slit, we 
neglected it. 

The combined frames were then flux-calibrated. To 
this aim, the photometric response of the instrument was 
measured with three standard stars, BD+28 4211, G191- 
B2B and GD71, observed with a 3.6" wide slit and through 
airmasses similar to the one of the observations of NGC 
588 (less than 1.2). We excluded the stellar lines from 
the measurement points, because of the difference of spec- 
tral resolution between our observations and the reference 
spectra. The atmospheric extinction function we adopted 
is the average one of the Observatorio de Roque de los 
Muchachos (La Palma, Spain), located at an altitude very 
similar to the one at CAHA. In each of the three spectral 
ranges, we found the response measures, as obtained with 
the different available spectra, to show the same chromatic 
trends, but constant discrepancies of ~0.1 magnitude. We 
also found oscillations of amplitude «3% on scales of 50- 
200 A, which we were unable to fit. Because of the drop 
in the instrumental response, the A < 4000 A range is 
more uncertain at ~15% with respect to the overall range 
B. Furthermore, in our data, a series of wide regions in 
the range Rl and, above all, R2 suffer severe extinction 
by telluric O2 and H2O unresolved molecular bands. Such 
extinction tended to cause numerical instabilities when 
fitting smooth functions on the measured photometric re- 
sponse, so we decided to correct them, even coarsely, to 
avoid these instabilities. For this, we synthesized the O2 
and H2O molecular bands, starting f rom individu a l lines 
as c atalogued in t he solar atlases of iMoore et alJ 1^966) 
and lMohlerl {l950). In both atlases, only wavelength cen- 
ters Xi and line strengths Wi are available, and we decided 
to model the bands as functions of magnitude loss of the 



form F(A;X) = N(X) £ s WiH{\ ~ A,), with X=0 2 or 
H2O, N(X), the relative column density of X on the line 
of sight, and H(X), the Gaussian spectral PSF of the spec- 
trograph. We then corrected the observed spectra of the 
standard stars for these extinction functions, using values 
of iV(Oa) and 7V(H20) that smoothed the corrected spec- 
tra best, and we computed the instrumental response in 
the ranges Rl and R2 with them. 

We checked the photometric consistency between the 
computed responses in the three ranges B, Rl, and R2. 
For this, we calibrated the three corresponding spectra of 
one of the standard source, BD+28 4211, observed in the 
same time interval, and searched for possible discontinu- 
ities from range to range. We found no such discontinuity, 
meaning that the three computed responses are mutually 
consistent from the point of view of absolute photometry. 

The last step in the calibration of the optical spec- 
tra was sky background removal. At each wavelength, this 
background was taken as a linear ramp fitted on two zones 
of the slit situated on both sides of the nebula and free 
of nebular or stellar emission. Figure [5] shows the optical 
spectrum along the whole wavelength range covered, inte- 
grated over the 81-pixel (45") zone of the slit where the 
extinction was computed. 

2.1.2. Estimation of the reddening and of the 
dereddened H/3 flux 

Once the spectra were reduced, we first performed an es- 
timation of the profiles of the color excess Eb-v(x) and 
of the dereddened H/3 flux Fp ol {x) along the slit, using 
the case-B line ratios of the hydrogen Bal mer decrement 
and a ssuming a gas temperature of 11000 K l(Vflchez et alJ 
1988j). We compared the fluxes of the Balmer lines Ha, 
H/3, H7, and H<5 at the different positions x along the slit 
common to both ranges B and Rl. The idea underlying 
the use of these four lines instead of two is the possibility 
of correcting the evaluations of Eb-v{%) and F^ eT (x) for 
the errors in the line fluxes. Such errors were split in three 
components: photometric residuals, underlying stellar ab- 
sorption lines, and data noise. 

We first measured the fluxes Fx(x) (X = a, /3, 7, S). 
For each line and at each position, the flux was measured 
as the sum of the pixels covered by the line, from which the 
continuum was previously removed. The fluxes were mea- 
sured on spectral windows wide enough to include the un- 
derlying stellar absorption lines entirely. Then, wc applied 
slight shifts to the line flux profiles, in order to account for 
object distortions in the frames (see Sect. l2.l"TJl . We also 
convolved the Ha profile by a mask so that its PSF along 
the slit matched the one of the three other lines. Finally, 
we dereddened the fluxes from foreground Galactic red- 
dening, using the Galactic extinction law llNandv et al 
| l975[ISeatoril979(l and E B -v = 0.045 l|Burstein fc Heiles 
Il984h . 

Using the LMC law l)Howarthlli983h for the extinc- 
tion due to M33, we then proceeded to correction of the 



4 



L. Jamet et al.: The giant Hn region NGC 588 



CM 

1 


1 20 


E 




o 


1 00 


in 




cn 




0) 


60 


m 




O 


40 


X 


20 






vL 





CN 


10 






u 


8 






( 71 


6 


o 






A 


o 






2 


X 








EH. 







4000 



5000 



6000 7000 
Wavelength (A) 



8000 



9000 



Fig. 2. Optical spectrum of NGC 588 shown in two flux scales. 



line fluxes for photometric residuals and for underlying 
stellar absorption lines, using the procedure explained in 
Appendix 1X1 Finally, from the corrected profiles F^? lr (x) 
(X — a, (3, 7), we proceeded to the fit of Eb-v{x) 
and Fp T {x). The E B -v{x) curve was smoothed with a 
Gaussian pattern of 6 pixels FWHM, because of signifi- 
cant noise in some zones of the slit. On the edges of the 
nebula, the smoothed curve was slightly negative so we 
clipped it to 0. Both E B -v(x) and F$ eT (x) curves are 
shown in Fig. |3| As an indication of the quality of the 
dereddening procedure, Fig. ^ shows the relative residu- 
als F™ rr (>)/FP rod (>) - 1, where F^ cd (x) is the flux of 
HX predicted from F^ er (x), Eb~v(x) and the theoretical 
HA/H/3 ratio, for the Ha, H/3, and H7 lines. Unexplained 
systematic residuals of a few % are found for H/3 and 
H7 around x = 250; but by changing artificially the line 
fluxes, we found that in this part of the slit, the uncer- 
tainty on Fp er (x) and Eb~v{ x ) is increased by less than 
2% and 0.01, respectively. 



measured. Each line was fitted by a Gaussian curve in a 
series of bins along the slit. The weakest lines were fitted in 
bins of several pixels whose sizes were ruled by their over- 
all signal-to- noise ratios (S/N), so as to ensure the conver- 
gence of the Gaussian curve fits. The fits were performed 
with the program TWODSPEC of the STARLINK 2 soft- 
ware suite; the [On] A3726+3729 doublet and of the [Ol] 
A6300 line, blended with an emission line of the sky, were 
measured with a personal IDL routine. Each bin was 
then dereddened using the corresponding value of Eb-v- 
Finally, to obtain the line intensities with respect to H/3, 
we summed the dereddened line fluxes in all the bins and 
divided the result by the dereddened H/3 flux summed 
over the same bins. In this procedure, we included the 
bins where Eb-v was clipped to 0. 



2.1.3. Line listing and measurement 



We established a list of detected lines from visual inspec- 
tion of the calibrated frames and selected the ones to be 



2 The STARLINK Project is run by CCLRC on behalf of 
PPARC. TWODSPEC was developed by T.N. Wilkins and 
D.J. Axon. 
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2.1.4. Detection limit and accuracy of the measured 
fluxes 

From measurements of the noise in the frames, we es- 
timated the 3cr detection limit of emission lines in the 
spectra. We found this limit to vary approximately from 
F(H/?)/2000 (A > 5000 A) to F(H/3)/500 (A » 3700 A). 
The S/N was then very high for the main lines. However, 
recombination lines of metals remained undetectable. In 
particular, the 3cr detection threshold for the Cn A4267 
and On A4651 was w F(H/3)/500. 

As already mentioned, in each spectral range, we found 
oscillations of ~3% amplitude in the residuals of the fit of 
the photometric instrumental response. These oscillations 
are even higher below A = 4000 A, causing a w 15% un- 
certainty in the photometric match of this range with the 
rest of the B range. Hence, when calculating line intensity 
ratios, we decided to add 3% of uncertainty in the fluxes 
of the lines with A > 4000 A, and 15% for lines with 
A < 4000 A. There are two exceptions to this rule: the H/3 
line (see Sect. I2.1.2|) 



and the ratios of lines belonging to 
a close doublet and then suffering identical photometric 
errors (e.g., [On] A3726/A3729, [Sn] A6731/A6716). 

Telluric absorption bands are another important 
source of uncertainty. The way we corrected the calibra- 



tion spectra for their effect is very approximate and let in 
an uncertainty of ^30% at A > 8200 A. However, the main 
source of error in this range comes from the lack of knowl- 
edge about the coincidence of the individual (unresolved) 
lines of the bands with the emission lines of the nebula, 
especially [Sin] A9069 and [Sin] A9532. Indeed, the flux 
of a nebular line can be severely reduced, depending on 
its exact position with respect to the telluric absorption 
lines. Consequently, we did not use the fluxes of the lines 
situated at A > 8200 A. 

Finally, an additional uncertainty in each line ratio 
results from a constant uncertainty of 0.02 in the Eg_y 
curve. Table|21presents the reddening-corrected intensities 
of all the lines measured in the spectrum and used in this 
work, in units of H/3 = 100, and the associated standard 
deviations (i.e. the uncertainties related to random errors 
only). 

2.2. Ground-based images 

We retrieved several narrow-band images of the ING 
archive. They were obtained through filters centered on 
Ha, H/3, [Om] A5007 and on the red side of the Ha con- 
tinuum, at the wavelength A = 6384 A. The Ha filter is 
also transparent to the [Nil] A65 (48+83) doublet, but the 
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Table 2. Dereddened fluxes of the lines and doublets used 
in this work, in units of H/3 = 100, and associated measure 
standard deviations in the same units. "Fluxes corrected 
for the SWS apertures. 



Line/Doublet 


Flux 


Standard deviation 


Optical slit 


[On] A372(6+9) 


141.54 


0.36 


[Nem] A3869 


36.86 


0.18 


[Sin] A40(69+76) 


1.58 


0.14 


[Om] A4363 


3.67 


0.07 


[Arrv] A4711 


0.26 


0.07 


[Arrv] A4741 


0.34 


0.05 


[Om] A5007 


404.79 


0.20 


[Nil] A5755 


0.20 


0.03 


Hei A5876 


11.94 


0.04 


[Sill] A6312 


1.38 


0.04 


[Nil] A6583 


10.69 


0.04 


[Sn] A6716 


12.73 


0.04 


[Sn] A6731 


9.17 


0.03 


[Arm] A7135 


9.47 


0.03 


[On] A73(20+30) 


3.09 


0.05 


Integrated fluxes 


[Om] A5007 


352.22 


12.46 


[Siv] All^m" 


37.81 


9.72 


[Nem] A16^m a 


56.51 


4.75 


[Sill] A19/im a 


81.42 


13.66 


[Sill] A33/xm a 


56.23 


9.73 
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214.47 


8.18 
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Fig. 5. Map of Eb-v established with the continuum-free 
Ha and H/3 images. The average uncertainty of this map is 
1=3 0.03. The computed reddening outside the ellipse cutout 
is less accurate, so it is not shown here. 



latter contributes little to the signal of the images and was 
neglected. 

Using some stars common to all the images, we sub- 
tracted the Ha continuum from the three other images, in 
order to keep only the signal of the nebular emission lines. 
We then established a reddening map, as shown in Fig. [5] 
This map was obtained after smoothing the Ha and H/3 
images with a 31x31 pixel mask, due to the small S/N 
ratio of the latter image. It can be seen that the redden- 
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Fig. 6. ISO-SWS detected lines. 



ing in the nebula is quite uniform, so we replaced it with 
a constant reddening. This simplification causes an error 
of only 1% in the computation of the total dereddened 
H/3 flux. We also established an [Om] A5007/Ha map, as 
shown in Fig. IB. 21 and discussed in Appendix iBl 

Knowing the position of the optical slit on the neb- 
ula and the H/3 flux measured throu gh it, and assum ing 
a distance of 800 kpc to the nebula ( Lee et a"fll2002l) . we 
calculated the total extinction-corrected H/3 luminosity of 
NGC 588: L p = 1.22 x 10 38 erg s" 1 . We tested this lumi- 
no sity by using the radio fl ux density at 1.4 GHz measured 
bv IViallefond et alJ {l986). The observed radio- to-H/3 ra- 
tio is F„(1.4 GHz) / Ff3 = ( 3 3±0.2) x 10~ 14 Hz" 1 W e used 
the e missivity formulae o f IPeauignot et al ] l|l99l|) for H/3 
and o f lOsterbrockl l|l974h for the radio continuum. In the 
latter, we included the free-free emissions induced by the 
H + and He + ions. Adopting a temperatur e of 11000 K 
and a n Hc + /H + abundance ratio of 0.083 fvilchez et all 
ljl98Sf) . we found a theoretical ratio F V (1A GHz)/^ = 
3.66 x 10~ 14 Hz -1 , a value higher by 10% than the ob- 
served one. The theoretical value of F V (1A GHz)/F ) g de- 
pends little on the temperature, and we can assess that 
the total dereddened H/3 flux that we found i s con sis- 
tent with the measurements of IViallefond et "all l|l986h to 
within 10%. Consequently, for comparison with the ISO 
data, we used our value of Lp with an associated uncer- 
tainty of 10%. 



2.3. ISO line spectra 

A number of infrared (IR) lines were observed with the 
Infrared Space Observatory (ISO), in the ranges of both 
the Large- and Small- Wavelength Spectrometers (LWS 
and SWS, respectively). We retreived 3 data sets from 
LWS observations and 1 from SWS observations (Table 

m 

The L WS data sets 8 80026 8 and 81601776 were an- 
alyzed by iHigdon et alJ l|2003f) . who found a value of 
3.51 x 10" 12 erg s" 1 cm" 2 for the [Om] A88^m flux. 
Given the important constraint provided by this line in 
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our study, we performed our own measurements on the 
three available data sets. For this, we used the ISAP 3 fa- 
cility. For each spectrometer and each detector (#4 and 
#5 for LWS01, #5 for LWS02), we obtained one spec- 
trum by selecting the scans to combine and by removing 
bad pixels. Even when processed, the LWS02 contained 
spurious features, so it was rejected. The measurements 
of the [Om] \88fim flux on the four LWS01 spectra gave 
compatible values. Combining them, we find a flux of 
(3.41 ±0.13) x 10~ 12 erg s -1 cm" 2 , in agreement with the 
value given bv iHiedon et al.1 ((2003) . In the following we 
adopted our value. Accounting for th e estimates given in 
the LWS Handbook l|Grv et al.ll2003|) on the repeatability 
of the LWS01 observations and the absolute photometric 
accuracy of this spectrometer, we estimate that the over- 
all uncertainty in the [Om] A88^m flux is 20%. Since the 
LWS aperture (84") is much larger than the nebula, no 
aperture correction was needed for the [Om] A88/im flux. 

The SWS spectra were processed with the ISAP fa- 
cility in the same way as the LWS data. We measured 
the fluxes of the lines [Neni] A16/mi, [Sin] A19^m, [Sin] 
A33/zm, and [Siv] All-in (Fig. 0). The SWS apertures 
range from 14" x 20" to 33" x 20" and are therefore smaller 
than the nebula. To estimate the total nebular fluxes in 
each of these lines, we assumed that their surface bright- 
ness distribution is identical to the one observed in the Ha 
image. This assumption is not entirely correct because of 
ionization stratification of the nebula. We tested it with a 
few photoionization models, using various ionizing spectra 
and density distributions, and found that it leads to an er- 
ror of w 10 % at most. Consid ering Table 5.3 of the SWS 
Handbook l|Leech et alJl2003|) . we adopted an instrumen- 
tal uncertainty of 20% in addition to the dispersions in 
the line fits and the error in the aperture corrections. 

The aperture-corrected intensities of the ISO detected 
lines are reported in Table [21 together with their non- 
instrumental uncertainties and the aperture correction 
factors. They are given relative to the total H/3 flux of 
the nebula (Sect. El}. 



2.0 2.5 
Log N, (cm- 3 ) 



Fig. 7. (T e ,log-/V e ) diagram for the following ra- 
tios. Temperature diagnostics (with vertical error 
bars): [Om] A4363/A5007 (full line, no symbol), 
[On] A73(20+30)/A372(6+9) (dotted line, square), [Sn] 
A40(69+76)/A67(16+31) (dot-dashed line, diamond), 
[Nil] A5755/A6583 (long-dashed line, circle), [Om] 
A5007/A88/im (short-dashed line, triangle). Density di- 
agnostics (with horizontal error bars): [Sn] A6731/A6716 
(full line, triangle), [On] A3726/A3729 (short-dashed line, 
square), [Sm] A19/mi/A34^m (dotted line, diamond). The 
error bars shown here are representative of the whole 
log N e /T e ranges. 

Table 3. Empirical temperatures derived from five line 
ratios 



Line ratio 



Temperature (K) 



[Om] A4363/A5007 11140 ± 180 

[On] AA7320+30/A3727 11210 ± 1180 

[Sn] AA4069+76/AA6716+31 9860 ± 940 

[Nil] A5755/A6583 11360 ± 840 

[Om] A5007/A88^m 8420 ± 690 



3. Empirical diagnostics 

3.1. Density and temperature 

A series of temperature and density diagnostics are avail- 
able from our observed line flux ratios. We established a 
(T e , log A e ) curve for five temperature and three density 
diagnostics, as seen in Fig.0 In this diagram, a 7% uncer- 
tainty was set to the [Sn] A6731/A6716 ratio, in order to 
account for the uncertainties in the collisional strengths of 
this ion. Likewise, the [On] A3726/A3729 ratio uncertainty 
was set to 15%. 

At temperatures close to 10000 K, a simple least- 
square fit of the [On] A3726/A3729, [Sm] A19/jm/A34^m 

3 The ISO Spectral Analysis Package (ISAP) is a joint de- 
velopment by the LWS and SWS Instrument Teams and Data 
Centers. Contributing institutes are CESR, IAS, IPAC, MPE, 
RAL and SRON. 



and [Sm] A19/im/A34/im line ratios with a single elec- 
tron density yields N e — 70 ± 60 cm -3 . However, while 
the [On] and [Sm] ratios indicate very similar densities, 
the [Sn] ratio indicates a smaller one. This is not nec- 
essarily inconsistent, since [Sn] lines are emitted further 
from the ionizing source than [Sm] and [On] lines, so that 
the electron density indicated by the [Sn] doublet is ex- 
pected to be smaller in case of an outward-decreasing gas 
density. Figure El shows the [Sn] A6731/A6716 and [On] 
A3726/A3729 ratio profiles along the slit with respect to 
the H/3 flux distribution. Both ratios are increasing func- 
tions of the electron density. While the [Oil] profile is 
nearly constant, the [Sn] one suggests a decrease in the 
outskirts of the nebula. However, in the low-density regime 
we observe here, interpretation of the [Sn] and [On] pro- 
files is not fully conclusive. The H/3 flux is another density 
indicator and we use it in Sect. Ifci . lTTl 
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Table 4. Empirical abundance ratios of various elements, 
with the ICFs specified 



20 30 40 50 60 70 
Relative position (pixels) 

Fig. 8. [Sn] A6731/A6716 and [On] A3726/A3729 ratio 
profiles along the slit. The H/3 flux profile is shown in 
the upper panel. 

Taking an electron density of 70 cm -3 , we derived em- 
pirical temperatures from the five corresponding diagnos- 
tics. The values are listed in Table [21 As seen, the four 
diagnostics using only optical lines, T opt (Oiii), T(On), 
T(Sii), and T(Nn), are in relatively good agreement with 
a mean temperature of « f 1000 K, T op t(Oni) being by far 
the most constraining of the four diagnostics. However, the 
[Oin] A5007/A88^m ratio gives a temperature Ti R (Oiii) 
that is lower by « 3000 K than T opt (Oiii). This discrep- 
ancy cannot be attributed to the error bars alone. Indeed, 
to yield a temperature of 11000 K, the [Om] A5007/A88/Ltm 
ratio would have to be 3.6 instead of 1.6. Such a large 
discrepancy cannot be attributed to errors of measure- 
ment and/or calibration. Instead, it suggests the presence 
of large temperature inhomogeneitics in the nebula. 

The discrepancy between the optical temperature di- 
agnostics and Tiri(Oin) might arise from two phenomena: 
(i) the optical slit samples a peculiar zone of plasma hotter 
than the overall gas of the nebula, making the tempera- 
ture diagnostics, representative of this zone higher than 



Ratio Value 


ICF 


Our work 


He/H (8.96 ± 0.42) x lO"^ 
O/H (1.45 ±0.27) x 10" 4 
N/O (4.12 ±0.95) x 10~ 2 
Ne/O (2.57 ± 0.14) x 10" 1 


He=He+±He ++ 

o=o+±o ++ 

N/0=N+/0 + 
Ne/0=Ne ++ /0 ++ 


Vi'lchcz ct al. (1988) 

r, l,r on.. 



O/H (2.0 ±0.3) x KT 4 
N/O (3.0 ±0.5) x 10" 2 
Ne/O (2.5 ±0.5) x 10" 1 



Tiri(Om) which is measured on the whole nebula; or (ii) 
temperature variations exist throughout the entire nebula. 
In the latter case, for each of the optical temperature di- 
agnostics, the more temperature-sensitive line (e.g., [Oin] 
A4363) significantly weights the hotter plasma zones, bi- 
asing the temperature measure toward a higher than av- 
erage value. This bias also affects Tir,(Oiii) but to a lower 
degree, hence making it smaller than the temperatures 
derived from optical line ratios. 

Assuming that the discrepancy between T opt (Oin) and 
Tir(Oiii) is not a sampling effect, we can estimate the 
average T and the mean square fluctuations T§ t 2 of 
the temperature in the ++ zone, both weighted by 
N e N(0 ++ ). By expanding to order 2 the sensivity of the 
lines to the temperature and assuming a density of 70 
cm" 3 , we found T = 8110±660 K and t 2 = 0.087±0.008, 
i.e. the temperature undergoes fluctuations of about 30% 
RMS. This amplitude of fluctuations is larger than those 
generally f ound in Hll regions (typical values of t 2 are 
0.02-0.05: lEsteban et alll2002|) . 

3.2. Abundances 

We computed a few empirical elemental abundance ratios, 
using N e — 70 cm -3 and the + and ++ optical tem- 
perature diagnostics. The ionic abundances were derived 
from lines measured in the optical spectra. The values and 
the ionization correction factors (ICFs) used are summa- 
rized in Tabled The error bars account for the uncertain- 
ties in the adopted density, in the empirical temperatures 
and in the line fluxes . The abundances ar e comparable to 
the ones inferred by IVflchez et al.l fmsft . We f ound the 
O/H v alue to be 30% of the solar one as given bv lLoddeTsl 
( 2003). At this stage, we did not compute the abundances 
of elements such as S or Ar, due to the uncertainties in 
their ICFs. 

The empirical elemental abundance ratios presented 
here served as a starting point for our modeling of the 
nebula. 
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Table 5. Observational values of the line ratios used to constrain the photoionization models and model results in the 
form (model-observation)/ (tolerance). The models are: [FS] Full Sphere; [HB] Hollow Bubble; [DDI] and [DD2] H/3 
constrained Density Distribution, 4.2 Myr spectrum; [DD2] H/3 constrained Density Distribution, 3.6 Myr Spectrum; 
[DF] Density fluctuations; [DG] Dust Grains. The geometric corrections for the models of Sect. and further sections 
are also written. a Values are valid for the whole nebula; the other data are representative of the optical slit. 
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[Om] A5007/H/3 
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-7.1 
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-5.9 


-0.8 


+0.9 


-1.5 


+2.2 


Nnl A6563/FOlll A(3726+91 


075±16% 


—0.1 


+0.0 


—0.4 


— 1.1 


— 1.0 


—0.6 


—0.3 


[Nein] A3869/[Om] A5007 


0.0911±5% 


-4.3 


-6.0 


-4.7 


-4.0 


-0.5 


-0.2 


+1.0 


[Arm] A7135/H/3 


0.0946±5% 


+1.7 


+1.7 
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+6.9 
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+5.6 
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-2.2 


-2.8 
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[Oi] A6300/H/3 


0.02±50% 
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[Om] A5007/H/T 
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-5.4 
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-4.6 


-0.9 


+0.8 


-1.1 
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Geometric corrections 
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1.07 


1.07 


1.40 


1.75 


1.35 


K 








1.48 


1.48 


0.52 


0.52 


0.12 


B 








0.58 


0.58 


1.45 


3.31 


1.15 



Stasifiska 



Stasinska 



4. General strategy for model fitting 

To m odel the nebula, we used the code PHOTO 
199CJ) with the update of atomic data listed in 
(120051) . The abundance ratios Mg/O, Si/O, S/O, Cl/O, 
Ar/O Fe/O a n d C/N were set to the solar values quoted 
from lLoddersI l 200.St. The input io nizing; spectra were 
those inferred by Jamet et alJ l(2004|) in their star-by-star 
analysis of the cluster, for an age ranging from 3.6 to 4.4 
Myr. 

The optical slit that we used is thin compared to the 
nebula, and we can consider that the line fluxes measured 
in the optical spectra are representative of a plane sec- 
tion in the nebula, perpendicular to the sky. Consequently, 
those lines were compared to model line intensities inte- 



grated in this plane slice. In the case of the line fluxes 
representative of the whole nebula ([Om] A5007/H/3 ratio 
from the images, lines observed with ISO), the model line 
intensities were integrated into the whole gas volume. 

Table summarizes the line ratios used to constrain 
the models, sorted according to the data they mainly di- 
agnose. The order of the lines in each ratio was chosen 
so that the ratio is an increasing function of the quantity 
it probes. We considered that satisfactory models should 
be able to reproduce all those constraints simultaneously. 
The observed values are quoted with their relative toler- 
ances in percentile of the observed values. For most of the 
line ratios, those tolerances were defined as the la measure 
uncertainties. The exceptions are the following. We set the 
tolerances on [Sir] A6731/A6716 and [On] A3726/A3729 to 
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7% and 15%, respectively, due to the uncertainties in the 
collisional strengths of the regarded transitions. We also 
increased by 25% the tolerance for the [Sm] A6312/[Sn] 
A67(17+31) and [Siv] All//m/[SlIi] A(19+34)^m diagnos- 
tics, in order to account for the poorly known dielectronic 
coefficients of sulfur. In Table the line ratios used as 
temperature diagnostics are followed by the correspond- 
ing empirical temperatures for a density N e — 70 cm~ 3 . 



5. Simple models 

5.1. Full sphere 

We first modeled the nebula in a very simple way. We 
assumed it to be a full sphere of constant density with 
a radius R = 1.5 x 10 20 cm which approximately fits the 
orientation-averaged size of the bright annulus that can be 
seen in the Ha image (Fig.^l. The input ionizing spectrum 
was the one obtained for an age of 4.2 Myr, the best-fit age 
from ljamet et all l|2004|) . We used the total H/3 luminosity 
of the nebula to constrain the relation between the filling 
factor e and the hydrogen density Nr. 

We considered several values of e. For e = 1, which 
leads to the smallest possible density, the model (labelled 
[FS] in Table EJ) produced [Om]/[Oii], [Sm]/[Sii] and 
[Siv]/[Siii] ratios that agree satisfactorily with the data 
but have too a small value of [Ariv]/[Anil]. For mod- 
els with e < 1, all four excitation diagnostic line ratios 
were underestimated. However, the main failure of the full 
sphere model is that it does not conciliate the temperature 
diagnostics based on optical lines and the [Om] optical/IR 
one. Indeed, it predicts a total amplitude of the variations 
of the electronic temperature T e of w 1000 K only. 

5.2. Hollow bubble 

The Ha image (Fig. shows that the nebula is better 
represented by a bubble than by a full sphere. Hence, we 
modeled the nebula as a hollow bubble whose thickness 
is 20% of its average radius R — 1.5 x 10 20 cm. We as- 
sumed e = 1, leading to a density = 11 cm -3 (de- 
rived from the H/3 luminosity) . As seen in Table this 
model, labelled [HB], fails to reproduce T opt (Oin), while 
it marginally agrees with Tir(Oiii). It also underestimates 
the excitation of the nebula. 

It seems clear that a more complex geometry for the 
gas is required in the models. Indeed, a complex density 
distribution is expected to change the ionization structure 
of the nebula, hence affecting the excitation diagnostics. 
Such a change would also modify the distribution of the 
main cooling ions (e.g., ++ ) and their efficiency to ra- 
diate energy away from the nebula. This may produce 
spatial temperature variations that are larger than in the 
case of a uniform density nebula. 



6. Models with constrained density distribution 

We used the well-observed distributions of the H/3, [On] 
A372(6+9) and [Om] A5007 flux profiles along the optical 
slit and the Ha and [Om] A5007 images to build a model 
of the gas density distribution. This model was used in all 
the photoionization models that follow. 

6.1. Construction of the density model 

6.1.1. RMS density assuming spherical symmetry 

In a nebula, the H/3 flux observed in a column of solid 
angle Q and coordinates (x, y) is: 

Q rZ 2 {x,y) 

Fp{x,y) = — / {ephvp N e N p e)(x,y,z)dz (1) 

47T J Zl {x,y) 

where z is the coordinate along the line of sight, ephvp, 
the H/3 emissivity, and N e and N p , the electron and proton 
density. Hence, assuming a given symmetry of the nebula, 
a model of the gas density distribution can be derived 
from the projected Fp(x,y) distribution. This is what we 
did for NGC 588, from the flux observed in the CAHA slit 
and from the Ha image. 

The Ha image suggests that the nebula is formed 
mainly by a bubble, possibly containing tenuous gas and 
surrounded by a halo. It also contains a bright knot lo- 
cated at the extremity of a filament. This knot and the 
other extremity of the filament are covered by the slit. 
With these considerations, we derived the following RMS 
density Nn y/e = f(x, z) in the plane covered by the opti- 
cal slit, from the long-slit profile Fp{x): 

— two "half-disks" located on either side of a central line 
of sight (of abscissa x = 0), each composed of a density 
plateau with f(x, z) — N% for r — \J x 2 + z 2 < R\, an 
annulus with f(x, z) = N2 > Nt for Rx < r < R2, and 
a fading halo with f(x, z) — N3 (R2/r) 2 for R2 < r < 
R3 

— two Gaussian round "blobs" on either side of the 
central line of sight, with the form f(x, z) = 
N G ex P (-((z - x G ) 2 + {z- z G ) 2 )/{2a 2 G )) ' 

This model is presented in Fig. EH where the good 
agreement between the observed Fp (x) profile and the one 
predicted from the density model can be appraised. In the 
model that we fitted, the values of Ni, N2, N3 and R3 
depend on the considered side of the central line of sight 
(x < or x > 0). We simplified this by adopting average 
values of those parameters in what follows; these values 
are listed in the 2nd column of Table We also neglected 
the smaller Gaussian blob, because of its small contribu- 
tion to the nebular flux (3% of the H/3 flux in the slit). 
Furthermore, we adopted a distance of 4 x 10 19 cm be- 
tween the bright blob and the stellar source, equal to the 
projected distance between this blob and the main con- 
densation of the stellar cluster. In what follows, we will 
call S the gas distribution that symmetrically surrounds 
the source, and G the bright Gaussian knot. G contributes 
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Fig. 9. Upper panel: H/3 long-slit profile (full line: ob- 
served; dashed line: model). Middle panel: model of the 
Nh y/e distribution (full line: gas surrounding the source; 
dashed lines: Gaussian blobs). In this panel, the value of 
N-r \fe for the two Gaussian blobs was divided by 4, for 
better visibility. Lower panel: schematic view of the model 
density profile in the plane section of the optical slit. The 
x axis is parallel to the slit, while the z axis is parallel 
to the line of sight. The two blobs (represented by black 
dots) may actually be located at any position along z. 



to 30% of the H/3 flux seen through the optical slit and 
cannot be neglected. 

We extended the model of density distribution to the 
case of the whole nebula, by using the Ha image (Fig.^J. 
For this, we stretched the spherical distribution S in the 
image plane, according to the elliptical shape of the bub- 
ble in the image, while conserving the original dimensions 



Table 6. Parameters of the model density distribution, 
fitted on the slit H/3 profile and on the image. The dis- 
tances are expressed in 10 20 cm, and the densities in 
cm~ 3 . "Projected minor and major semi-axes. ^Projected 
distance to the source. 



Parameter Slit value Image value 



Ri 


1.01 


0.97x1.61" 


Rz 


1.24 


1.19x1.98° 


R 3 


2.49 


2. 39x3.99° 


Ni 


4.81 


6.13 


N-> 


9.13 


10.9 


N 3 


5.09 


6.04 


Rg" 


0.40 


0.40 




0.17 


0.17 


N G 


40.3 


40.3 



along the line of sight. We then re-fitted N\, N2 and N3 
on the parts of the Ha image unoccupied by the bright 
central filament. The values obtained, listed in Table 
are slightly different from the ones obtained for the opti- 
cal slit. G contributes to « 7% of the total H/3 flux of the 
nebula. 

The knot G is relatively confined and is not centered 
on the source. Since our code only processes spherically 
symmetric gas distributions, we modeled G with a shell of 
Gaussian density profile, but we integrated the line fluxes 
in a way adapted to its non-radial gas distribution. 

6.1.2. Accounting for non-sphericity 

The model presented above assumes spherical symmetry 
of its different components, in contradiction with what 
stands out of the Ha image. Indeed, the projected im- 
ages of both S and G are elliptical, showing that their gas 
distributions are not isotropic. Furthermore, such asym- 
metry can exist along the line of sight with respect to the 
plane of the sky and it that case, it remains undetectable 
because of the projection of the signal received from the 
nebula. 

Our photoionization code assumes spherical symmetry, 
meaning that the accurate geometry of S cannot be pro- 
cessed. This is why we decided to model S with a spherical 
gas distribution obtained by homothctic transform of the 
three-component distribution TVh e = f(r) established in 
Sect. l6.1Jl This transform, illustrated in Fig. ^3 is meant 
to infer sizes to the different structures (e.g., the bubble) 
that are averages of the real dimensions of these struc- 
tures. We used two homothetic factors, A for the lines seen 
through the optical slit and B for the lines integrated over 
the whole nebula. A and B are basically different, because 
they are representative of different portions of the nebula. 

In Sect. IbM.ll we assumed G to be spherical. 
Nevertheless, it may be stretched or squeezed along the 
line of sight with respect to the plane of the sky. In this 
case, the volume of G is not the one that we inferred in 
Sect. 16.1.11 Since its H/3 flux is fixed, its density is also 
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Fig. 10. Illustration of the geometrical correction for the 
gas distribution of S. The arrows indicate the direction of 
the observer, (a) Model inferred in Sect. l6.l!Tl (b) Possible 
shape of the real distribution, (c) Approximation of (b) 
adopted in photoionization models in Sect.EH3 

different from the one initially inferred. This effect must 
be taken into account, because of the importance of the 
density on the ionization state of the gas. Hence, we in- 
troduced a factor K of stretching G along the line of sight 
with respect to the projected view. 

The geometric corrections that we introduced have sig- 
nificant consequences for the ionization state predicted 
by the models. In particular, they modify the [Om] 
A5007/[On] A372(6+9) and [Om] A5007/H/3 ratios in both 
S and G. We used those ratios to constrain the coeffi- 
cients A, B and K , following the procedure explained in 
Appendix iBl 

6.2. Model results 

Model [DDI] of Table[S]was established with a filling fac- 
tor of e = 1. All the parameters other than the density are 
the same as those of model [FS] and [HB] (Sect.EJ. While 
Tra(OlIl) is marginally reproduced, T op t(0ni) is still much 
underestimated. 

In order to attempt to fit the optical temperature di- 
agnostics — and temporarily ignoring the optical/IR di- 
agnostic — we modified the input carbon and sulfur abun- 
dances, since these two elements are the most efficient 
coolants among those whose abundances were not directly 
measured. The other parameters, including A, K and B, 
were unchanged. We were able to fit the temperature di- 
agnostics only with very small carbon and sulfur abun- 
dances: even if the latter are set to 1/10 of the original 
values, T O pt(0m) and T(Sn) are almost too small com- 



pared to the observed ones (column [DD2] of Table |3J). 
Such small abundances are physically unlikely, so we re- 
jected the model. 

7. Further refining of the models 

7.1. Using the degree of freedom on the stellar ionizing 
spectrum 

So far, we only used the 4.2 Myr cluster model of 
Ijamet et al. (2004j), which corresponds to their best fit 
model. However, acceptable cluster ages range from 3.6 to 
4.4 Myr. The hardness of the Lyman continuum decreases 
with increasing age. Therefore, we used the cluster spec- 
trum corresponding to 3.6 Myr to test whether it is at least 
able to explain the 11000 K temperature diagnosed by the 
optical lines (ignoring again the optical/IR diagnostic). 

We discarded the filling factor e = 1, because it re- 
quires values of A and B that are too large, which is equiv- 
alent to assuming that the nebula is extremely stretched 
along the line of sight. For e = 0.1 and using the original 
carbon and sulfur abundances, T op t(Oiii) is still underes- 
timated. However, this issue is resolved if the sulfur abun- 
dance is divided by 2 (model [DDS]), a value that can be 
considered acceptable. 

Although using an ionizing spectrum harder than the 
one first used improves the prediction of the optical tem- 
perature diagnostics, the large temperature discrepancy 
between the latter and Tir,(Oiii) remains unexplained. 

7.2. Appending local density fluctuations 

Most GHRs present complex local structures, such as fila- 
ments and co ndensations (e.g., NGC 604 a nd 30 Dor: see, 
respe ctively, iMaiz-Apellaniz et alJ 120041: IWalborn et alJ 
l2002f) . Such structures may exist in NGC 588 and remain 
invisible in our ground-based images due to insufficient 
angular resolution. Local structures in the gas are usually 
modeled by using filling factors e < 1, i.e. assuming that 
the gas is entirely condensed in uniform clumps occupying 
a fraction e of the available volume. However, it is more 
realistic to assume that the density varies continuously, 
being higher in the bright structures than in the regions 
situated between them. 

A continuous density distribution with non-zero den- 
sity between the bright structures may increase the ampli- 
tude of the spatial temperature variations since, due to a 
different distribution of the cooling ions, the temperature 
in the diffuse component is expected to be different from 
that in the filaments. 

7.2.1. Formalism of the fluctuations 

The RMS density of Sect.Elwas multiplied by the following 
periodic radial function: 



s{r;H,r ) = 



L+(l-L)|cos(7rr/r )r 



H 



(2) 
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Fig. 11. Radial profiles of iVn, T e , and the relative abun- 
dances of the most coolant ions, in model [DF] . 



where L = min(s) / max(s) defines a diffuse gas component 
between a series of clumps, whose condensation level is 
described by the parameter n; H is a normalization factor 
such that < f 2 > r = 1. H ensures that the observed H/3 
flux distribution will not be modified within scales larger 
than r since, locally, the H/3 is proportional to s(r) 2 . We 
note that max(s 2 ) = 1/H. Furthermore, H only depends 
on L and n, and the actual value ro is only required to 
be small, in order to simulate small-scale gas structures. 
Finally, we will deal not with L, but with <fi = L 2 /H, the 
fractional contribution of the diffuse component to the H/3 
flux. We applied the density fluctuations only to S, and 
kept the original smooth density profile for G. 

7.2.2. Results 

We computed a few models with various values of H and 
0, using the stellar spectrum for 3.6 Myr and the standard 
elemental abundances. The different models gave similar 
results. The one shown in Table [5] ([DF]) was obtained 
for H = 0.1 and </> = 0.2. It can be seen that the fluctu- 
ations of density do not resolve the discrepancy between 
T opt (OlIl) and Tir(Oiii). 



Fig. 12. Radial profiles of Na, U, T e (three upper panels), 
and fractional [Oiii] A5007 and [Oiii]A88/im fluxes inte- 
grated through increasing circular apertures (lower panel) , 
in model [DG] . 



Figure ^2 shows the radial distribution of N^, T e) 
0+/0, 0++/0 and H°/H in S. The variations caused 
by the density fluctuations to the ionization state induce 
variations in the cooling rate of the different chemical el- 
ements. In particular, in the outer parts of the nebula, 
neutral hydrogen becomes an efficient cooling agent by 
collisional excitation of Lya. Unfortunately, the resulting 
variations of temperature between the density clumps and 
the surrounding gaps are smaller than 500 K. 

Close to the source, the temperature fluctuations asso- 
ciated to the clumps are due in large part to the ionization 
structure of sulfur. Consequently, we multiplied the sulfur 
abundance by 2 and checked the changes in the temper- 
ature distribution. We found that this distribution was 
modified very little, so we concluded that density fluctua- 
tions are not a solution to the T op t(Oiii)-T IR (Oiii) conflict. 

7.3. Models with dust 

The effects of dust in the thermal balance of plan etary 
nebulae were studied bv lStasifiska fc Szczerbal ( 200ll) with 
detailed models. They show that dust grains provide ad- 
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ditional heating by photoelectric emission, proportional 
to the dust-to- hydrogen mass density ratio Pd/pn and to 
the local ionization parameter U(r) = Q(H°)/(47rr 2 N e c). 
Hence, in a nebula where the dust density and the aver- 
age ionization parameter are large enough, the decrease in 
U (r) with distance r to the source causes a significant neg- 
ative gradient of the temperature as a function of r. Local 
density fluctuations also change the temperature signif- 
icantly toward lower values in clumps due to a smaller 
value of U(r) in them. 



We adopted the grain size distribution of lMathis et alJ 
l)l977|) for half the total grain mass and a small grain pop- 
ulation (with grain sizes < 0.01 /im) for the other half 
mass. We limited our study of dust grain effects to the 
case of the cluster spectrum for 3.6 Myr, which gives the 
highest value of Q(H°) among our set of cluster spectra. 
This allowed us to increase U(r) as much as possible. It 
also allowed us to maximize the threshold of Pd/pn above 
which dust grains would absorb too large a fraction of 
the ionizing photons and cause the ionization bound to be 
closer to the source than the observed bound. This thresh- 
old is Pd/pn = 5 x 10~ 4 for the adopted ionizing stellar 
spectrum. 

We adopted H = 0.01 and <fi — 0.5 for the density 
fluctuation function. With these parameters, the diffuse 
gas and the clumps contribute equally to the H/3 flux and 
U(r), which is inversely proportional to s(r), changes by 
a factor as large as l/^/tpH — 14 between the maxima of 
the clumps and the diffuse gas component. 

Column [DG] of Table [S] shows the line ratios pre- 
dicted with pdj pa = 5 x 10~ 4 . Even with such a high 
dust-to-gas density ratio, we were unable to obtain tem- 
perature fluctuations large enough to explain the conflict 
between T opt (Oni) and Tir(Oiii) and both temperatures 
are largely overestimated. Indeed, in addition to a global 
increase of temperature, the three effects of dust on the 
thermal balance are a significant steepening of the tem- 
perature in a small zone close to the source, which only 
contributes to w 10% of the [Om] A5007 and [Om] A88/zm 
fluxes, a drop of the temperature in the faint outskirts of 
the nebula and a drop of the temperature in the clumps 
by up to w 2000 K, which is insufficient to explain the 
temperature discrepancy between the observed diagnos- 
tics. This is illustrated in Fig. where the radial pro- 
files of Ns., T e , and U(r) are plotted together with the 
fractional fluxes of the [Om] A5007 and [Om] A88/^m in- 
tegrated through circular apertures whose radii are given 
by the abscissa axis. 

To reproduce at least T op t(Om), Pd/pa would have to 
be smaller than in model [DG]. In that case, by being 
smaller than in the latter, the amplitude of the temper- 
ature variations would be even further from accounting 
for the T op t(Oni)-TiR,(Oin) discrepancy. Dust grains are 
therefore an unsatisfactory solution. 



20 



10 



o 
o 
Ifi 
u 



1 



-20 




-20 -10 10 
Aa (arcsec) 



20 



Fig. 13. Zones of extraction of spectra in Sect. 17.41 



7.4. Spatial distribution of the stars 

It is well known that in a specific point of a nebula, 
the energy gain provided by photoionization mainly de- 
pends of the ha rdness of the local Lyman continuum (e.g., 
IStasihskall200^ . and in particular on the temperature T* 
given by: 



3 S Vn °°0--voMFua v du 



L n Fv<J v j{hv)dv 



(3) 



where F v is the flux received from the stars, Hvq the hydro- 
gen ionization threshold, and o v the hydrogen ionization 
cross-section. In the case of NGC 588, the ionizing stars 
are distributed on spatial scales comparable to the size 
of the nebula, and one may expect the Lyman continuum 
received to vary from point to point, as it is harder close 
to hotter stars. In that case, the temperature structure of 
the gas is sensitive to the distribution of the stars. 

To evaluate the impact of the heterogeneity of the clus- 
ter on the temperature structure of the gas, we divided 
the cluster into three zones, indicated in Fig. 1131 and ex- 
tracted the stellar spectrum of each of them, adopting the 
age of 3.6 Myr for the cluster. As shown in Tabled the 
three zones individually produce a significant fraction of 
the ionizing photons and have different temperatures T*. 
Altogether they provide all the ionizing flux of the cluster. 
The zones SE, close to the knot G, and SW, also covered 
by the optical slit, are dominated by the two Wolf-Rayet 
st ars (WR) of the clus ter, i.e. respectively stars #2 and #1 
of lJamet et all ll2004|) . They are the brightest and hottest 
two stars of the cluster. Zone N is dominated by a group 
of three stars, much cooler and less luminous than the two 
WR. 

We re-computed the model [DDS] of Sect. 17. ll but re- 
adopted the original sulfur abundance, and successively 
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Table 7. some spectral properties of the three zones N, 
SW and SE and of the whole cluster, assuming an age of 
3.6 Myr. Qtot(H°) is the total photon luminosity of the 
cluster in the Lyman continuum range. 



Zone Q(H u )/Q tot (H u ) T, (kK) 

N 25% 34 

SW 26% 38 

SE 49% 48 

All 100% 41 



Table 8. Empirical temperatures yielded by the spectra 
of the zones N, SW and SE and of the whole cluster. 



zone 


Topt(Oiii) 


T(Oii) 


T(Sn) 


T(Nn) 


Tir(Oiii) 


N 


9990 


9860 


7960 


10090 


9580 


SW 


10220 


10080 


8070 


10300 


9760 


SE 


11610 


11450 


8830 


11610 


10880 


All 


10840 


10740 


8430 


10920 


10250 



replaced the spectrum of the cluster by the ones of the 
three zones N, SW, and SE, each rescaled to the ionizing 
rate of the whole cluster. The results are shown in Table 
[8] in terms of empirical temperatures at N e = 70 cm~ 3 . 
While the spectra of N and SW result in very similar tem- 
peratures of the gas, the one of SE yields temperatures 
higher than the latter by ss 1500 K. Thus, close to the 
stars of the zone SE, the gas temperature may undergo a 
significant excess with respect to the overall nebula. This 
is the case of G, in particular, much closer to SE than to 
the other parts of the cluster, at least in the projected sight 
of the nebula. Hence, we may a priori suspect G of being 
responsible for an increase in T op t(Oni), since the latter is 
measured through the optical slit, where G provides ~40% 
of the [Oin] A5007 flux. However, by replacing the model 
of G inferred from the total cluster spectrum with the one 
derived from the spectrum of WE, we oberved an increase 
of T O pt(0iii) by 300 K only. 

The spatial structure of the cluster of NGC 588 turns 
out not to provide a satisfactory explanation of the 
r op t(OlIl)-Tm,(OlIl) conflict that we encountered. 

8. Sources of energy other than the ionizing flux 

8.1. Shocks 

A hypothesis that is often advocated to explain the large 
temperature inhomogeneities observed in some GHRs is 
the supply of energy in the form of shocks resulting from 
the winds of massive stars and from supernovae. Those 
shocks are expected to raise the temperature abruptly in 
thin slabs of gas. However, alteration of the thermal struc- 
ture induced by shocks in nebulae has been studied lit- 
tle (s ee, however, the case of NGC 2363: iLuridiana et alJ 
2001). A "ho t-spot" model of this pheno menon has been 
proposed by iBinette fc Luridianal l)2000|) , which we ap- 
plied it to NGC 588. 



No supernova remnant was identified wit h in NG C 588 
cither in the X-ray sur vey of IPietsch et al. | ll2004) or m 
the radio-optical one of lGordon et al.l (ll QQflh . We there- 
fore assumed that the kinetic power released by the stellar 
cluster is conveyed only by winds. Using the theoretical 
wind powers implemented in Starburst99 l)Leitherer et alJ 
Il999l) . we found that the stars of NGC 588 release a ki- 
netic power E = 1.1 x 10 38 erg s _1 , close to the H/3 lu- 
minosity of the nebula. On the other hand, the photoion- 
ization gain of the nebula is 15-20 times the H/3 lumi- 
nosity. As a result, if we assume that all the energy of 
the shocks enhances the c ollisionally excited lines of the 
nebula, the term Thoat of IBinette fc Luridianal l)2000|) is 
0.05-007. Given the low metallicity of the nebula, this 
leads to temperature fluctuations below t 2 = 0.01. This is 
clearly insufficient to resolve the T op t(Oiii)-TiR(Oiii) con- 
flict. To reach t 2 — 0.087, the kinetic power of the stellar 
winds would have to be at least 5 times larger than that 
computed. 

8.2. Electron conduction from hotter plasma layers 

Hot plasma slabs generated by shocks in a nebula are 
expected to convey heat to the surrounding gas by 
electron conduction, creating transition zones between 
those slabs and the gas excited by photoionization alone. 
iMacieiewski et all ( 19961) constructed a model nebula con- 
sisting of gas sheets, each of width H, separated by thin 
layers of very hot and tenuous gas. They studied the ther- 
mal structure of the sheets undergoing heat conduction 
by electrons, and showed that this structure is able to en- 
hance the [Oni] A4363 line, while lines such as H/3 and 
[Oin] A5007 are modified little. Through their Eqs. (9), 
(10), and (14)-(18), they suggest a relation between H, 
the density N p , the temperature T oq in zones undergo- 
ing photoionization alone, and the temperature diagnostic 
Topt(Oni). Adopting T oq = 8000 K and ignoring Eq. (18), 
we found H ~ 10 15 -10 16 cm, within the range of densi- 
ties N p inferred in Sect. 16.1.11 This estimate is of course 
rough, because it does not accurately take into account the 
density structure of the nebula (smooth or with clumps), 
the ionization state of the gas, and the real value of T eq . 
However, given the size of the nebula, our estimate of TL 
holds the presence of ~ 10 4 -10 5 layers of very hot plasma; 
even with higher values of H., the number of those layers 
would be much too large to be realistic. We can therefore 
reject the electron conduction of heat as an explanation 
of the large difference between T op t(Oni) and Tir(Oiii). 

9. Uncertainties in O/H and Ne/O 

Since the thermal balance of NGC 588 is not fully under- 
stood, the estimates of abundance ratios, either empirical 
or obtained with models, are uncertain not only through 
the line ratios that are used to probe them, but also 
through the discrepancy that exists between the different 
temperature diagnostics. Wc estimated this uncertainty 
for the ++ /H and Ne ++ /0 ++ abundance ratios, as ++ 
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Table 9. ++ /H + and Ne ++ /0 ++ abundance ratios derived from different diagnostics. The error bars include the 
uncertainties related to the temperature diagnostics and to the [Om]/H/3 and [Neiii]/[Oin] line ratios that we used. 



Set of data Temperature Line ratio Abundance ratio 
0++/H+ 

Optical slit r op t(Oiii) = 1140 ± 180 K [Om]A5007/H/3 = 4.05 ± 0.18 (1.04 ± 0.07) x 10~ 4 

Integrated r m (Om) = 8420 ± 690 K [Oiii]A5007/H/3 = 3.52 ± 0.20 (2.55 ± 0.72) x 10~ 4 

Ne ++ /0 ++ = 

Optical slit r opt (Oiii) = 1140 ± 180 K [Neni] A3869/[Om]A5007 = 0.0911 ± 0.0047 0.257 ± 0.014 

Integrated Tir(Oiii) = 8420 ± 690 K [Nem] A16/LHn/[Om]A88/Ltm= 0.257 ± 0.085 0.164 ±0.059 



and Ne ++ are the two ions for which we have observations 
in both optical and IR ranges. We compared the empirical 
++ /H and Ne ++ /0 ++ ratios obtained by using the long- 
slit values of [Om] A5007/H/3, [Nem] A3869/[0m] A5007, 
and T opt (Oiii) with those derived from the integrated val- 
ues of [Om] A5007/fb3, [Nem] A16/xm/[Om] A88^m, and 
Tir(Oiii). 

The results are reported in Table El It can be seen 
that the computed ++ /H abundance ratio depends a lot 
on the selected temperature diagnostic and that the ac- 
tual value of ++ /H is known to within ±0.2 dex. Since 
++ is the dominant ion of O in the nebula, the elemen- 
tal abundance ratio O/H is uncertain by a similar amount. 
Assuming that the empirical ++ /0 derived from the op- 
tical slit data is correct (0 ++ /0=1.39), and depending 
on the adopted [Om] temperature, O/H ranges approxi- 
mately from 1.5 x 1CP 4 to 3.6 x 10~ 4 . The two estimates 
of the Ne ++ /0 ++ ratio, to which Ne/O is expected to be 
close, differ by a smaller factor of 1.6. 

It is interesting to note that the value of Nc ++ /0 ++ 
derived from Tir(Oiii) and [Nem] A16/zm/[Om] A88/xm 
closely matche s the solar Ne/O ratio recommended by 
lLoddersl lf2003). which is 0.15. This agrees with the usual 
idea that Ne/ O is close to s olar in environments such as 
GHRs (e.g., iGarnetd Eoolfr . Furthermore, the O/H es- 
timate derived from Tir(Oiii) and the integrated [Om] 
A5007/H/3 ratio is found to be close to the average one of 
supergiant stars of M33 sit uated at similar galact ocentric 
distances (O/Hra 3 x 10" 4 : llJrbaneia et"aHl2005|) . Hence, 
the integrated optical/IR line ratios suggest "standard" 
abundances in the nebula. On the contrary, the optical 
spectrum suggests peculiar abundances with a ~50% O/H 
depletion and a Ne/O excess of ~70 %. 

10. Conclusion 

This work was motivated by the frequent use of giant Hn 
regions to determine the chemical composition of the inter- 
stellar medium in galaxies. The usual methods are simple 
in principle but need to be validated by realistic models. 

We have carried out a comprehensive study of the giant 
Hn region NGC 588. We gathered and carefully combined 
a large number of observations consisting of spectroscopic 
data in the UV, visible, and infra-red domains, as well as 
of narrow-band images. The most outstanding property of 



this Hn region is that the temperature derived from the 
[Om]A5007/A88/im ratio is lower by about 3000 K than 
the temperature obtained with four different optical tem- 
perature diagnostics. This suggests either the presence of 
large temperature inhomogcncities or the incorrectness of 
flux measurements. After a careful analysis of the obser- 
vational data, we excluded this latter possibility, because 
it would imply that the [Om] A88/J,m line flux has been 
overestimated by a factor of 2. We then constructed pho- 
toionization models with the aim of reproducing the obser- 
vational data, including all the temperatu re diagnostics. 

In a previous paper ijjamet et al.ll2004j) , we conducted 
a detailed star-by-star analysis of the cluster responsible 
for the ionization of the Hn region and obtained an ion- 
izing spectral energy distribution based on the most up- 
to-date stellar atmosphere models. We used this spectral 
energy distribution as input to photoionization models tai- 
lored to reproduce all the observed spectral diagnostics. 
Our procedure takes the different apertures into account 
through which the spectroscopic data are obtained. We 
started with the simplest density structures, i.e. a homo- 
geneous sphere and then a bubble. Indeed, most of our 
understanding of giant Hn regions is based on such mod- 
els. Since these simple models failed to reproduce all the 
temperature diagnostics, we proceeded to more complex 
representations. Our strategy was to stick as closely as 
possible to the observational constraints provided by the 
spectra and images and to explore the free parameters 
left systematically. These include abundances of uncon- 
strained elements such as carbon, density fluctuations, 
dust grains, spatial extent of the ionizing cluster, shock 
heating and conductive heating. No satisfactory solution 
was found. We stress that, although there are other cases 
of Hn regions for which published photoionization models 
were not able to reproduce the observed temperatures, this 
is the first time that such a thorough analysis of a giant 
Hn region has been performed and such a large number 
of observational constraints used. 

The fact that no acceptable model has been found yet 
has important consequences. Two explanations can be ad- 
vanced: 

— Either the energy balance of this object is not un- 
derstood. This t hen might also apply to giant Hn re- 
gions in general ijEsteban et alJl2002t and references 
therein). In that case, abundances derived from opti- 
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cal lines are not as accurate as usually thought, even in 
objects where temperature diagnostics are available. 
— Or, the uncertainty in the [Oni] A88/im fluxes mea- 
sured by ISO is larger than expected. If this is the 
case, results based on the analysis of this line, even in 
other contexts, should be regarded with caution. 

For the time being, in the case of NGC 588, we can only 
say that the O/H ratio probably lies somewhere between 
w 1.5 x 10~ 4 and 3.6 x 10~ 4 . However, we found that the 
O/H and Ne/O ratios derived from optical/IR integrated 
data are more consistent with at least some comparison 
sources (supergiants in M33 for O/H, the Sun for Ne/O) 
than those derived from the optical spectrum and should 
probably be preferred to the latter. 

Alternatively to temperature fluctuations, abundance 
inhomogeneities might be a solution to the conflict be- 
tween the temperature diagnostics. They have been advo- 
cated in planetary nebulae to explain the inconsistencies 
between abundances derived fr om optica l forbidden lines 
and from recombination lines Such conflicts 

are a lso observed in giant Hn regions (e.g., lEsteban et alJ 
l2002h . albeit to a lesser extent. Though abundance inho- 
mogenei ties are theoretically expected not to occur in Hn 
regions llTeno rio-Taele 1996), the current view i s changing 
( Stasihska et al.l2005l:lTsamis & Peauignotl2005() . A mea- 
surement of oxygen recombination line intensities in NGC 
588 would provide an important constraint to investigat- 
ing the presence and significance of abundance inhomo- 
geneities in this object. In the case where there the latter 
is present, one needs to investigate the precise meaning of 
the chemical abundances derived with various techniques. 
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£,x{x) = ± ax{x). The equivalent widths of the stellar 
lines are approximatel y the same for Ha, H/3, and H7 (e.g., 
lAuer fc Mihalasiri972l) . but may vary spatially if the stel- 
lar content is heterogeneous. Consequently, we described 
the stellar contribution to the flux errors as the sum of 
the contributions of different subsets i of the cluster, each 
being associated to an equivalent width Wi and continua 
F^ x (x). The sum of the errors 6x(x) is: 



Sx(x) 
F x (x) 



F x {x) F x {x) 



(A.4) 



The linear equation IA.3I applies if the relative errors 
6x(x)/ Fx{x) are small, i.e. if the photometric errors ex 
are much smaller than 1, if the equivalent widths of the 
nebular lines with respect to the stellar continuum are 
much greater than the ones of the stellar lines, and if the 
signal-to- noise ratios (S/N) Fx{x)/ax{x) are high. 



Appendix A: Correction of Balmer lines used for 
dereddening 

Let us consider a set of hydrogen Balmer lines, whose 
fluxes are measured at a series of positions x along a slit. 
At a given position x, the measured value of the flux of 
"HA", F x (x), differs from its real value F^ l (x) by a 
measure error 8x(x) = Fx(x) — F x 2al (x). Furthermore, 
F x al (x) is related to Fg eal (x) by the following formula: 



2.5 1ogFf al (z) 



2.5 log(^ cal (x)i? x ) 
+E B _v(x) Up- fx) 



(A.l) 



where Rx is the theoretical HA/H/3 intensity ratio, and 
f x the extinction law. By applying this equation to Ha, we 
can derive the following formula, especially for the bright 
H7 and H<5 lines: 



(fx-U) In-Fr'M 

+(/„-/*) In (A.2) 

Analogously, we can write the flux of HA predicted 
from the observed fluxes of Ha and H/3, namely F x lcd (x). 
Due to the measure errors, F x ed (x) and Fx(x) are dis- 
crepant, and if the errors are small compared to the fluxes, 
then the discrepancy term Ax(x) = Fx(x) — F x cd (x) is 
given by: 



{fp - fa) 



Ax(x)-S x (x) 
F x (x) 



{fa - fx) 



Hfx - U) 



5p{x) 
Ffs(x) 
6 a (x) 



F a (x) 



(A.3) 



We split the flux errors into three terms: spatially in- 
variant residuals ex of the photometric calibration, stellar 
absorption features underlying the nebular Balmer lines, 
and the usual line flux dispersions (such as photon noise) 



A.l. Photometric errors 

Let us consider the set of positions x unaffected by stellar 
absorption lines, and where the nebular lines were mea- 
sured with high S/N. Combining Eas. IA.31 and IA.41 give: 



{fp ~ U) 



with 



Axjx) 
F x (x) 



{fx - fp)e a + {fa - fx)ep 

+ (//3 - fa)tX ± T X {x) 



(A.5) 



■_{x) = {fx~fp? 



v a {x) 
F a {x) 



+ {U - fx? 



t s2 j' (J x{x) 



<?p{x) 
F {x). 

(A.6) 



The contribution r\x of the photometric residuals to 
Ax{x) I Fx{x) can be estimated with the following for- 
mula: 

J: x {Ax{x)/Fx{x))/t x (x) 



Vx 



(A.7) 



(A.6 



•qx relates ex to e a and ep: 

fx — fp fa — fx 
£x — Vx - -? T^a 7 T e P- 

fp - J a fp - J a 

Since they arc residuals of the fit of the instrumental pho- 
tometric response, we can write ex = Oi&x- Furthermore, 
if they are not correlated between them, we can assume 
they minimize the error balance \ 2 — Ex e x Ivx 1 m ore 
specifically given by: 



! (e Q , £/3 ) - 



x 



1 ( fx - fp ^ fa — fx ^ 
4>X W/3 - fa 6a fp ~ fa 



'VX 



(A.9) 



The estimations and uncertainties of e a , ep, and more gen- 
erally ex can be easily derived from this linear x 2 an d from 
Eq. IA.8I The correction of the line fluxes for these errors 
consists in multiplying them by 1 — ex- 
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A. 2. Underlying stellar absorption lines 

Making use of fluxes corrected for photometric errors and 
measured with high signal-no-noise ratio, we can write, for 
H 7 : 



ft f ^^zM 



Wiu>i(x) ± t 7 



= (/ 7 - 9p)^44 + if" - fi) 



F a (x) 
+(fe ~ /a) 



F 7 (x) 



(A.10) 
(A.H) 



The estimate of the equivalent widths Wi is given by the 
minimization of the following linear % 2 : 



x 2 = 



Appendix B: Setting the geometric corrections A, 
B, and K 

We used the well-observed long-slit profiles of H/3, [Oiii] 
A5007, and [On] A372(6+9), and the [0m] A5007/Ha map 
to constrain the coefficients A, B and K. From those data 
(Fig. IB.l1 and lB~2|) . it is evident that the knot G is more 
excited than S, even though its RMS density is higher 
than the latter. This can be easily explained by its short 
distance to the stellar source (Sect. 

Coefficients A and B strongly influence the predicted 
ionization state of S, while K acts on the ionization state 
of G. We decided to set them, for each model, with the 
following method. 

We manually separated the contributions of S and G 
to the long-slit fluxes of [Om] A5007 and [On] A372(6+9), 
and derived the values of [Om] A5007/[Oil] A372(6+9) for 
both components. We found F([Om])/F([Oii]) = 2.14 ± 
0.15 for S and 5.68 ± 0.65 for G. 

For each model, we first computed the predicted fluxes 
of the lines observed through the optical slit. Unless oth- 
erwise mentioned, we set the coefficients A and K to 
reproduce the [Om] A5007/[On] A372(6+9) ratios of S 
and G. Once those coefficients were set, we reported the 
ratio R between the predicted and the observed value 
of [Om] A5007/H/3. Then, we calculated the model line 
fluxes integrated over the whole nebula. In the latter cal- 
culation, coefficient B was chosen so as to reproduce the 
model/observation ratio R, but for the integrated lines 
this time. 

The geometric transforms performed on the gas den- 
sity distributions (/(r) for S, g(r) for G) were accompa- 
nied with a re-normalization of the latter to conserve the 
predicted H/3 fluxes. Indeed, the integrated volume of gas 
in S is proportional to A 2 for the optical slit and to B 3 for 
the whole nebula, while the volume of G is proportional 
to K. Since the local H/3 emissivity is proportional to the 
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Fig. B.l. Profiles of H/3 (upper panel), [Om] A5007/H/3 
(middle panel), and [Om] A5007/[On] A372(6+9) (lower 
panel) along the optical slit. 



square density, the transform of f(r) consisted in replac- 
ing it by f(r/A)/A (optical slit) or f(r/B)/B 3 / 2 (whole 
volume), and g(r) was replaced by g(r) / K 1 / 2 . 
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Fig.B.2. Map of the nebular [Om] A5007/Ha ratio. The 
ratio is too noisy to be shown outside the ellipse cutout. 



